%% Basal Traction paper 2012 Apr
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\RequirePackage[fleqn]{amsmath}
\mathindent0pt

\documentclass[tc]{copernicus}

%% Language %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\usepackage[english]{babel} %francais, polish, spanish, ...
%\usepackage[T1]{fontenc}
%\usepackage[ansinew]{inputenc}
%\usepackage{lmodern} %Type1-font for non-english texts and characters
%
%%% Packages for Graphics & Figures %%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{graphicx} %%For loading graphic files
\usepackage{float} %%For loading graphic files
\usepackage{subfig} %%Subfigures inside a figure
\usepackage{boldtensors}

%
%%% Math Packages %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\usepackage{amsmath}
%\usepackage{amsfonts}
%\usepackage{amssymb}
%\usepackage{amsthm}
%
%%% Line Spacing %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%\usepackage{setspace}
%%\singlespacing        %% 1-spacing (default)
%%\onehalfspacing       %% 1,5-spacing
%%\doublespacing        %% 2-spacing
%
%
%%% Other Packages %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\usepackage{a4wide} %%Smaller margins = more text per page.
%\usepackage{fancyhdr} %%Fancy headings
%\usepackage{longtable} %%For tables, that exceed one page

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% DOCUMENT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}

\title{The effect of errors in bed elevation on adjoint based assimilation of
surface velocity}


\author[1]{Kevin Sack}
\author[2]{Jesse V. Johnson}


\affil[1]{Department of Mechanical Engineering, CERECAM University of Cape Town }
\affil[2]{University of Montana Department of Computer Science}

%% The [] brackets identify the author to the corresponding affiliation, 1, 2, 3, etc. should be inserted.

\runningtitle{~TEXT}
\runningauthor{~TEXT}

\correspondence{Name\\ (EMAIL)}

%\received{}
%\pubdiscuss{} %% only important for two-stage journals
%\revised{}
%\accepted{}
%\published{}

%% These dates will be inserted by the Publication Production Office during the typesetting process.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\firstpage{1}
\linenumbers
\maketitle
\begin{abstract}
A finite element numerical model is used to explore the interaction between the
combined basal geometry and traction with the surface velocity for a section of
a terrestrial ice sheet. The model domain is a two-dimensional flowline profile
extending far from the ice divide outwards. The model compares results between a
reference basal geometry and a series of significantly perturbed counter parts.
Finally an optimization routine is employed on the basal traction field, to
compensate for the perturbations and to minimize the difference between
perturbed and reference surface velocities.  
\end{abstract}


\introduction
Over recent years, the study of glaciers has received increased attention in many academic pursuits. Through environmental concerns, glaciers pose a potential risk to us over geological time scales; through rising sea levels via large scale melt water contributions or alternatively through ice advancement on bordering civilizations. In Climatology, glaciers hold, trapped beneath layers of preserved ice, the key to understanding historical atmospheric conditions (including chemical composites, global temperatures and flow patterns). Finally, as a global resource, glaciers hold $60\%$ of the world's fresh water. The understanding of glaciers is an important pursuit best achieved through numerical models.\\


Glaciers exist in areas where there is a net accumulation in snow on a yearly basis. As each year brings another layer of snow the previous layers are compacted through gravitational forces. Snow undergoes a transformation into solid ice over a period of several years. Ice is usually considered as a solid but is more accurately defined as a highly viscous, non-linear fluid, requiring a certain yield stress before undergoing flow \citep{pgm}.\\

\subsection{Paper focus}
A problem that effects glacier models is the inaccuracy of available data to verify modelling results. In particular interest to this paper is the collection of bedrock elevation data and characterising the bed geometry. The measurement of bedrock elevation is taken with specific radar sounders developed for ice. The quality of the data is limited by:
begin{itemize}
	\item The nature of radar technology (ie the echo interpreted as bedrock may not originate at the point directly under the sensor).
	\item The spacing between measurement track lines.
	\item The method employed to assimilate data between the measured track lines.
\end{itemize}


Comparatively, data collected for ice surface velocity can be captured more accurately and at higher resolution, which sets surface velocity as one of the benchmark measures of glacier model accuracy. This project is motivated by the hypothesis that given high resolution surface velocity data, we can compensate the influence of inaccurate bedrock elevation data, or bed geometry, by manipulating the basal friction, through the application of a scalar field defined on the bed geometry.\\


This would be set up in the form of an optimization problem, where basal friction would be changed to minimize the error between a high quality ``known'' ice surface velocity and perturbed counter parts.

%\subsubsection{Presentations of paper}
%
%\begin{itemize}
%	\item Chapter 2 will outline the properties of ice, introduce the governing equations and describe the non-linear rheology of ice flow.
%	\item Chapter 3 will cover the introduction of boundary conditions, the weak formulation of the problem and the Finite Element numerical method employed.
%	\item Chapter 4 will outline the terrain generation method employed to create a realistic bed geometry. 
%	\item Chapter 5 will show some basic results and initial analysis on our reference model.
%	\item Chapter 6 will introduce the optimization routine used and the mathematical formulation of the optimization problem.
%	\item Chapter 7 will show the optimization results and analysis.
%	\item Chapter 8 will cover the concluding remarks and describe the future directions for this model.
%\end{itemize}

\section{Ice Sheet Modelling}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Properties of Ice}
Ice can crystallize in nine different ways under various combinations of pressure and temperature. The temperature and pressure combination in nature only allows for the crystalline structure \textit{Ih}. Each oxygen atom bonds with four additional oxygen atoms, with two hydrogen atoms in between each bond. $H_{2}0$ molecules attach to each other in a hexagonal mesh, that forms a \textit{basal plane} and each of these planes lie on top of each other. The higher number of bonds within the basal plane compared to the bonds \textit{between} the planes means that when a stress is applied, the molecular structure favours deformation like a deck of cards. To a lesser extent gliding is also possible in the prismatic and pyramidal planes. These basic deformations are illustrated below:
 
 \begin{figure}[H]
  \centering
  	\subfloat[Basal]{\includegraphics[width=0.24\columnwidth]{pictures/Mol_def1}}                
  	\subfloat[Prismatic]{{   }\includegraphics[width=0.25\columnwidth]{pictures/Mol_def2}}
  	\subfloat[Pyramidal]{{   }\includegraphics[width=0.25\columnwidth]{pictures/Mol_def3}}\\
  \caption{Glide planes in the Hexagonal ice-Ih crystal}
  \label{fig:glide}
\end{figure}


\subsection{Governing Equations}
Ice behaves like a very slow moving liquid, with flow laws governed by Navier-Stokes equations and the incompressibility condition. The high viscosity, $\eta$ and low velocities $~v$ characterise ice flow as ``creeping flow'', allowing for the Navier-Stokes equations to be reduced to the simpler Stokes equations. The dependency of $\eta$ on the strain rate makes ice flow highly non-linear.

\subsubsection{Incompressibility}
The structure of solid ice is practically incompressible. The assumption of the incompressible condition is a good approximation on the model, considering the negligible upper layer of compacting snow. The incompressibility equation is given by:

\begin{equation}
\textbf{div}(~v) = 0\\
\label{eq:continuity}
\end{equation}

Where $~v$ is the velocity vector, describing the motion of the ice. 

\subsubsection{Conservation of Momentum}
Conservation of momentum, for a fluid, is obtained by an application of Newton's second law and is given by:

\begin{equation}
\rho \left( \frac{\partial ~v}{\partial t} + ~v \nabla \cdot ~v \right) = \textbf{div}(~\sigma) + ~f\\
\label{eq:momentum}
\end{equation}

Where $\sigma$ is the Cauchy stress tensor. The Navier-Stokes equations are obtained from (\ref{eq:momentum}) by assuming that the Cauchy stress, $~\sigma$, is the sum of a constitutive viscous term (proportional to the gradient of velocity) and pressure. In our case of semi-static slow moving flow, the inertial terms can be neglected, while the assumptions on $~\sigma$ are upheld, resulting in what is commonly called Stoke's flow:

\begin{equation}
~0 = \textbf{div} (~\sigma) + ~f\\
\label{eq:stokes}
\end{equation}

Where the force is defined by gravitational acceleration $g = 9.81 m.s^{-2}$ downwards and the density of ice $Ih$, $\rho=911 kg.m^{-3}$:
\begin{equation}
 ~f = -\rho ~g
\label{eq:force}
\end{equation}

\subsection{Stress and Strain}
In three dimensions, using the Cartesian coordinates, we have our velocity vector defined as, $~v = \{u, v, w\}$. The Cauchy stress tensor $~\sigma$ and the strain tensor $~\varepsilon$ are defined as:\\
\begin{equation}
\begin{array}{cc}
	~\sigma = \left[ \begin{array}{ccc}
		 \sigma_{xx} & \sigma_{xy} & \sigma_{xz}\\
	 	 \sigma_{yx} & \sigma_{yy} & \sigma_{yz}\\
		 \sigma_{zx} & \sigma_{zy} & \sigma_{zz} \end{array} \right] ,
  &  ~\varepsilon = \left[ \begin{array}{ccc}
		 \varepsilon_{xx} & \varepsilon_{xy} & \varepsilon_{xz}\\
	 	 \varepsilon_{yx} & \varepsilon_{yy} & \varepsilon_{yz}\\
		 \varepsilon_{zx} & \varepsilon_{zy} & \varepsilon_{zz} \end{array} \right] 
\end{array}
\label{eq:stress}
\end{equation}

Note that $~\sigma$ is symmetric. We will often be dealing with the strain rate tensor, also know as the rate of deformation tensor, $\dot{~\varepsilon}$, which is given by (and related to) the velocity gradients: 


\begin{equation}
\dot{~\varepsilon} = \frac{1}{2}\left( \nabla ~v + \left( \nabla ~v \right)^T \right)
\label{eq:straindot}
\end{equation} 

Following from the incompressibility condition (\ref{eq:continuity}) it is clear that

\begin{equation}
\textbf{tr}(\dot{~\varepsilon}) = \dot{\varepsilon}_{xx} + \dot{\varepsilon}_{yy} + \dot{\varepsilon}_{zz} = 0
\label{eq:compress}
\end{equation}

\subsubsection{Yield Stress}
In some materials there are no deformations at stresses below a certain stress, called a \textit{yield stress}. In other materials, such as ice, deformations are so small for low stresses that a theoretical yield stress is assumed. For ice sheet modelling, the driving stress that deformations are dependant on result from hydrostatic pressure. To be more specific the key dependency for deformation to occur require a significant pressure difference inside the glacier to drive flow. Horizontal flow will only occur if there is a slope to the ice surface that will induce a pressure difference within the material.\\


When we model glacial flow, instead of creating a yield stress condition, we ensure a significant surface slop is prescribed, in this way we avoid the complexity of modelling low stress behaviour. We impose the following restriction for surface slope $\alpha$:

\begin{equation}
\alpha \neq 0
\label{eq:surface}
\end{equation} 

\begin{figure}[H]
	\centering
		\includegraphics[width=0.6\columnwidth]{pictures/pressure.png}
	\caption{The ice deforms under its own weight, a gradient is needed to induce an internal stress difference between points that would drive horizontal flow}
	\label{fig:pressure}
\end{figure}

\subsubsection{Deviatoric Stress} 
The stress tensor can be broken down into the volumetric and deviatoric parts, which is essential to the formation of Navier-Stokes and Stokes equations. The deviatoric component described a material's deformation and as such will be of more interest to us. The deviatoric component, $\acute{~\sigma}$, is described by:

\begin{equation}
\acute{~\sigma} = ~\sigma - p\textbf{I}
\label{eq:dev}
\end{equation} 

Where $p$ is the mean normal (or spherical) stress:

\begin{equation}
\nonumber p = \frac{1}{3}\left(\sigma_{xx} + \sigma_{yy} + \sigma_{zz} \right)
\end{equation}

\subsection{Constitutive relation}
The most predominantly used constitutive law in glaciology, is a relation between stress and strain through a viscosity parameter:

\begin{equation}
\acute{\sigma}_{ij} = 2\eta \dot{\varepsilon}_{ij}
\label{eq:constit}
\end{equation}

To further develop the relationship between stress, strain and flow, we need to introduce effective stress and effective strain rate $\sigma_{e}$ and $\dot{\varepsilon}_{e}$ respectively, the measures introduced below are derived from the square root of \textbf{\textsl{II}}, the second scalar invariants of the respective tensors. 

\begin{eqnarray}
\sigma_{e} &=& \frac{1}{\sqrt{2}}\left( \acute{~\sigma} : \acute{~\sigma}\right)^{1/2} \label{seff}\\
\dot{\varepsilon}_{e} &=& \frac{1}{\sqrt{2}}\left( \dot{~\varepsilon} : \dot{~\varepsilon}  \right)^{1/2} \label{steff}\\
\nonumber
\end{eqnarray}


Using the invariants to describe the constitutive relation ensures that ice flow is independent of the orientation of the coordinates. The above definitions also correspond to the von Mises criteria of yielding in the theory of plasticity. Note that the incompressibility condition was used to obtain (\ref{steff}). Furthermore, we introduce a flow law, that $\dot{\varepsilon}_{e}$ is a function of  $\sigma_{e}$ \citep{Nye1}:

\begin{equation}
\dot{\varepsilon}_{e} = f(\sigma_{e})
\label{eq:flow}
\end{equation}

Note, since $\sigma_e$ is defined by terms of the deviatoric stress, $\acute{~\sigma}$, this definition constructed above has excluded the direct effect of hydrostatic pressure on flow. For a homogeneous, isotropic material, the above described flow law applies to all points at all times. In a more general case, that is beyond the scope of this paper, the above flow law can also be made dependant on temperature, moisture content and ice fabric.\\

Now, using (\ref{eq:constit}), (\ref{seff}) and (\ref{steff}) one can write:

\begin{equation}
\sigma_{e}^{2} = 4 \eta^{2} \dot{\varepsilon}_{e}^{2}
\label{eq:flow1}
\end{equation}
\begin{equation}
\Rightarrow \sigma_{e} = 2 \eta \dot{\varepsilon}_{e}
\label{eq:flow2}
\end{equation}

Taking into account (\ref{eq:flow}),
\begin{equation}
\Rightarrow \eta = \frac{\sigma_e}{2f(\sigma_e)}
\label{eq:flow3}
\end{equation}

Which finally allows us to write:
\begin{equation}
\acute{\sigma}_{ij} = \frac{\sigma_e}{f(\sigma_e)}\dot{\varepsilon}_{ij}
\label{eq:flow4}
\end{equation}

We approximate the flow function by a power law relation, \citep{Glen}:

\begin{equation}
f = \left( \frac{\sigma_e}{B} \right)^n
\label{eq:}
\end{equation}

Which properly allows us to express the constitutive relation, commonly know as Glen's flow Law as:
\begin{equation}
\acute{\sigma}_{ij} = \frac{\sigma_e^{1-n}}{B^n}\dot{\varepsilon}_{ij}
\label{eq:glen1}
\end{equation}

This constitutive flow law additionally relates deviatoric stress to velocity gradients (\ref{eq:straindot}). When combined with our conservation of momentum equation it allows us to fully express the Navier stokes and stokes equations.\\

We will often write Glen's flow law in the following form which is easier to manipulate:
\begin{equation}
\dot{\varepsilon}_{ij} = A \sigma_e^{n-1}\sigma_{ij}\\ \ \ \footnote{In the numerical scheme $\sigma_{e}$ has a tiny regularization parameter added to it, to prevent a singularity at zero shear rates}
\label{eq:glen2}
\end{equation}

Where $A = (1/B)^n$ and is know as the fluidity parameter. For \textit{Ih} ice moderately low temperatures ($\approx 263K$), the regularization parameter $n$ has been measured to be $n=3$ and $A$ has been measured to be $1 \times 10^{-16}Pa^{-3}yr^{-1}$.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Formulation}

A 2D representation of ice sheet flow, represents movement along a flow line without any loss of generality, given that a glacier's width does not change and there is no convergence from other flows into our ice sheet. In our case, a 2D form is particularly desirable when studying the impact the bed geometry has on surface velocity. It allows for clear analysis, which at a later stage could be generalized to the 3D representation. Perhaps the largest advantage of the 2D analysis is the computational efficiency which allows for a large number of model simulations.\\

 The following assumptions have been taken in our model of ice sheet flow:\\
\begin{itemize}
	\item Constant temperature has been assumed throughout the ice sheet.
	\item The ice body is taken to have uniform density and behave like an isotropic and continuous material.
	\item The effects of basal melting are ignored in this model.
	\item Friction is assumed to be uniform at the ice-bedrock interface for the initial reference set up. In the later chapters of the project, we rescale friction as the optimization variable.
	\item The bedrock does not deform under stress and none of the ice is absorbed into it.
	\item The ice flow is not dynamic. We limit the experiments to the steady state case.
\end{itemize}

\begin{figure}[H]
\centering
\includegraphics[width=.75\columnwidth]{pictures/setup}%
\caption{2D representation of ice sheet flow}%
\label{fig:setup}%
\end{figure}

Formulation of the problem will follow Voigt notation. The following are vector quantities:

\begin{equation}
\begin{array}{ccccc}
   ~\sigma = \left[ \begin{array}{c}
		 \sigma_{x} \\
		 \sigma_{y} \\
		 \sigma_{xy}  \end{array} \right] ,
  
  &  \nabla_s = \left[ \begin{array}{cc}
		 \frac{\partial}{\partial x} & 0 \\
		 0 & \frac{\partial}{\partial y}\\
		\frac{\partial}{\partial y} & \frac{\partial}{\partial x} \end{array} \right] ,	

  &	~v = \left[ \begin{array}{c}
		 v_{x} \\
		 v_{y}  \end{array} \right] ,\\
		 
  &  \dot{~\varepsilon }= \nabla_s ~v = \left[ \begin{array}{c}
		 \dot{\varepsilon}_{x} \\
		 \dot{\varepsilon}_{y} \\
		 2 \dot{\varepsilon}_{xy}  \end{array} \right] ,

  &	~I_s = \left[ \begin{array}{c}
		 1 \\
		 1 \\
		 0  \end{array} \right]		 		 
		 
		 \end{array}
\label{eq:Voigt}
\end{equation}


\subsection{Boundary Conditions}
For our specified domain $\Omega$, Dirchlet and Neuman Boundaries follow the notation $\Gamma_{D}$ and $\Gamma_{N}$ respectively. We have three distinct regions on the boundary that need to be constrained: The surface, the bed and the vertical boundary at each end of our 2D ice sheet.\\

The surface maintains a stress free boundary condition:
\begin{equation}
~\sigma \cdot ~n = ~0
\label{eq:surface}
\end{equation}

Where $~n$ is the outward normal unit vector. For the boundary condition at the bed, we take a simplified sliding law, similar to that originally defined by Weertman \citep{fowler}. We make the constraint that the velocity at the bed remains tangential to the boundary:

\begin{equation}
~v \cdot ~n = 0
\label{eq:noupslide}
\end{equation}

Furthermore, the Shear stress at the bed is related to the tangential velocity by the frictional forces acting at the rock-ice interface: 

\begin{equation}
~\sigma \cdot ~n = \beta^2 ~v 
\label{eq:bedslide}
\end{equation}

Where $\beta^2$ is a positive scalar representing the magnitude of frictional forces at the bed and due to (\ref{eq:noupslide}), the velocity along the bed moves in the direction, $\hat{~n}$, the unit tangent to the bed. This enables us to resolve the values for $~\sigma \cdot ~n$ along $\hat{~n}$. \\

To maintain our assumption of conservation of mass, we impose periodic boundary conditions on the velocity $~v$ at either end of our Ice sheet model. The following precautions have been set up in the model to ensure the periodic boundary conditions will not interfere with our experiments:

\begin{enumerate}
	\item Ice thickness at each end is identical
	\item The bedrock slope and geometry is kept uniform at each end 
\end{enumerate}

\begin{figure}[H]
\centering
\includegraphics[width=.75\columnwidth]{pictures/setupPBC}%
\caption{2D representation of Ice sheet flow, with periodic boundary conditions}%
\label{fig:setup2}%
\end{figure}

The above boundary conditions can be condensed into the following expressions:
\begin{align}
&~v &=\bar{~g}  &\;\;\parbox{5cm}{on $\Gamma_{D} $} \label{bc1}  \\
&~\sigma \cdot ~n &= \bar{~h} &\;\; \parbox{5cm}{on $\Gamma_{N} $} \label{bc2} \\
\nonumber
\end{align}

\subsection{Weak formulation }

Starting with the strong form of the problem:

\begin{eqnarray}
-\nabla^{T}_s ~\sigma &=& \bar{~f} \label{s1}\\
\textbf{div}(~v) &=& 0 \label{s2}\\
\nonumber
\end{eqnarray}

In the weak formulation, we cast the strong form into an equivalent form that is suited for solving with the Finite Element Method. This is achieved by by obtaining the inner product of the governing equations with a set of weighting, or test, functions. The solution to the weak formulation is, under certain conditions, equivalent to the solution to the strong form \citep{hughes}.\\

We will be using a mixed velocity-pressure formulation to achieve this. Given the set of trial velocity solutions $V$ and test functions $W$ which are defined as:

\begin{eqnarray}
U &=& \left\{ ~v | ~v\in \left[ H^1(\Omega) \right]^2,~v = \bar{~g} \text{ on $\Gamma_{D}$} \right\} \\
W &=& \left\{ ~w | ~w\in \left[ H^1(\Omega) \right]^2,~w = 0 \text{ on $\Gamma_{D}$} \right\} \\
\nonumber
\end{eqnarray}

Additionally we introduce pressure test functions $~q$ and trial functions $~p$ from the set $P$.

\begin{equation}
P = \left\{ ~q, ~p| ~q,~p\in\left[ L^2(\Omega) \right]^2\right\} \\
\nonumber
\end{equation}

The weak form for (\ref{s1}) is given simply by pre multiplying the equation by an arbitrary pressure test function $~q$ and integrating it:

\begin{align}
\int_{\Omega}{~q^T \textbf{div}(~v)} \textnormal{d}\Omega &= 0 \\
\label{Weak2}
\end{align}

Given $\bar{~f}$, $\bar{~h}$ and $\bar{~g}$ we can manipulate the strong form by pre multiplying (\ref{s1}) by an arbitrary test function $~w$ and integrating it. This yields:

\begin{align}
\nonumber -\int_{\Omega}{~w^T \nabla_s ~\sigma} \textnormal{d}\Omega &= \int_{\Omega}{~w^T\bar{~f}} \textnormal{d}\Omega \\
\nonumber
\end{align}

Applying Green's theorem to the first term:

\begin{align}
\nonumber \int_{\Omega}{ (\nabla_s ~w)^T ~\sigma} \textnormal{d}\Omega &= \int_{\Omega}{~w^T\bar{~f}} \textnormal{d}\Omega + \oint_{\Gamma}{~w^T ~\sigma \cdot ~n} \textnormal{d}\Gamma \\
\Rightarrow \int_{\Omega}{ (\nabla_s ~w)^T ~\sigma} \textnormal{d}\Omega &= \int_{\Omega}{~w^T\bar{~f}} \textnormal{d}\Omega + \oint_{\Gamma_N}{~w^T \bar{~h}} \textnormal{d}\Gamma \label{Weak1}\\
\nonumber
\end{align}

In the last line the boundary conditions and restriction to the test function $~w$ on $\Gamma_D$ was enforced. We know recount our definition of deviatoric stress (\ref{eq:dev}), which can be expressed in Voigt notation as:

\begin{equation}
~\sigma = \acute{~\sigma} + p ~I_s \\
\label{eq:devVoight}
\end{equation}

Applying our Glen's flow law (\ref{eq:glen1}), allows us to write this as:

\begin{align}
~\sigma = \frac{\sigma_e^{1-n}}{B^n} \left[ \begin{array}{ccc}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1/2 \\
\end{array} \right] \dot{~\varepsilon} +p\left[ \begin{array}{c}
1 \\
1 \\
0 \\
\end{array} \right] 
\label{eq:ConVoigt}
\end{align}

In a compact, more convenient form this will be written as:

\begin{align}
\nonumber ~\sigma &= \bar{~D} \dot{~\varepsilon} + p ~I_s\\
~\sigma &= \bar{~D} \nabla_s ~v + p ~I_s
\label{eq:ConVoigt2}\\
\nonumber
\end{align}

It is important to note that $\bar{~D}$ is a function of $~v$, characterising this as a non-linear problem. Substituting (\ref{eq:ConVoigt2}) into (\ref{Weak1}), gives us the following:

\begin{align}
\int_{\Omega}{ (\nabla_s ~w)^T \bar{~D} \nabla_s ~v} \textnormal{d}\Omega + \int_{\Omega}{ (\nabla_s ~w)^T p ~I_s} \textnormal{d}\Omega \ ... \nonumber \\ = \int_{\Omega}{~w^T\bar{~f}} \textnormal{d}\Omega + \oint_{\Gamma_N}{~w^T \bar{~h}} \textnormal{d}\Gamma \label{Weak2}\\
\nonumber
\end{align}

Noticing that $(\nabla_s ~w)^T p ~I_s = \textbf{div}(~w)^T p$ allows us to write this as:

\begin{align}
\int_{\Omega}{ (\nabla_s ~w)^T \bar{~D} \nabla_s ~v} \textnormal{d}\Omega + \int_{\Omega}{ \textbf{div}(~w)^T p} \textnormal{d}\Omega  \ ... \nonumber \\ = \int_{\Omega}{~w^T\bar{~f}} \textnormal{d}\Omega + \oint_{\Gamma_N}{~w^T \bar{~h}} \textnormal{d}\Gamma \label{Weak2}
\end{align}

Which allows us to write a single combined statement for the mixed velocity-pressure weak form:

\begin{align}
\int_{\Omega}{ (\nabla_s ~w)^T \bar{~D} \nabla_s ~v} \textnormal{d}\Omega +\int_{\Omega}{~q^T \textbf{div}(~v)} \textnormal{d}\Omega \ ...\nonumber \\ + \int_{\Omega}{ \textbf{div}(~w)^T p } \textnormal{d}\Omega =   \int_{\Omega}{~w^T\bar{~f}} \textnormal{d}\Omega + \oint_{\Gamma_N}{~w^T \bar{~h}} \textnormal{d}\Gamma \label{Weak}
\end{align}

Our aim is to to find $~v(x,y) \in V$ and that satisfies (\ref{Weak}) for all $~w \in W$ and for all $p \in P$.\\


\subsection{Finite Element formulation}
The Finite Element Method approximates solutions to the governing equations based upon their weak formulations. We solve the system by discretizing the domain into a mesh composed of elements and nodes, as illustrated below in figure (\ref{fig:femmesh}) \footnote{Image adapted from \citep{femmesh1}}. For this project triangular elements were chosen. The global approximations, spanning the entire domain will be referred to with the superscript $h$ (e.g $~v^h$ for velocity).\\
 
 \begin{figure}[H]
  \centering
  	\subfloat[Domain of interest]{\includegraphics[width=0.45\columnwidth]{pictures/femmesh1}}                  	   						\subfloat[Discretized domain]{{   }\includegraphics[width=0.43\columnwidth]{pictures/femmesh2}}\\
  \caption{Illustration of Finite Element discretization of a T shaped beam. Elements are represented by lines and nodes correspond to all intersections of elements.}
  \label{fig:femmesh}
  \end{figure}


On a specific element $e$ we will use the notation $~v^e$, and assume that $~v^e$ is non-zero only on element $e$. We express the approximation of each element by interpolating the nodal values of each element by the use of Shape functions (also known as basis or interpolation functions). An elemental approximation for $~v^e$ with $n_{en}$ elemental nodes can be written as:

\begin{equation}
~v^e = \sum^{n_{en}}_{i}{N^{e}_{i} \varphi^{e}_{i}}
\label{eq:element}
\end{equation}

For convenience, we will compose this in a matrix form and follow the basic method of \citep{fish}:

\begin{align}
~v^e = ~N^{e} {~\varphi}^{e}
\label{eq:element2}
\end{align}

These shape functions fulfil the ``kronecker delta''property; each Shape function $N_i =1$ at node $i$ and zero at other nodes. In this way, the nodal values are interpolated throughout the element \footnote{Interpolation can be linear, quadratic or higher order depending on the chosen order of the Shape functions}. In a similar manner, $~p^e$ can also be approximated, using the collection of Shape functions $\tilde{~N}^e$.\\

\begin{figure}[H]%
\centering
\includegraphics[width=0.7\columnwidth]{pictures/Shape}%
\caption{Quadratic Shape functions acting over a 1D three node element}%
\label{shape}%
\end{figure}

The global approximations are obtained by gathering the contributions from each element. For a mesh of $n_{el}$ elements, we can write: 
\begin{align}
\nonumber ~v 	& \approx ~v^h = \sum^{n_{el}}_{e}{~N^{e} {~\varphi}^{e}} \\
\nonumber ~w 		& \approx ~w^h = \sum^{n_{el}}_{e}{{~N^{e}~c^{e}}} \\
\nonumber ~p 		& \approx ~p^h = \sum^{n_{el}}_{e}{\tilde{~N}^{e}~p^{e}} \\
 		~q 		& \approx ~q^h = \sum^{n_{el}}_{e}{\tilde{~N}^{e}~q^{e}}
\label{eq:approx}
\end{align}

Further more, 

\begin{align}
\nonumber \nabla_s ~v^h &= \sum^{n_{el}}_{e}{\nabla_s ~N^{e} ~\varphi^{e}} = \sum^{n_{el}}_{e}{~B^{e} ~\varphi^{e}} \\
\nonumber \nabla_s ~w^h &= \sum^{n_{el}}_{e}{\nabla_s ~N^{e} ~c^{e}} = \sum^{n_{el}}_{e}{~B^{e} ~c^{e}} \\
 \textbf{div}(~w^h) &= \sum^{n_{el}}_{e}{~I^{T}_{s}~B^{e} ~c^{e}} = \sum^{n_{el}}_{e}{~L^{e} ~c^{e}} 
\label{eq:approx2}
\end{align}

Taking the weak form as a system of equations, grouped by like constants:
\begin{align}
\int_{\Omega}{ (\nabla_s ~w)^T\bar{~D}\nabla_s ~v} \textnormal{d}\Omega + \int_{\Omega}{\textbf{div}(~w)^T ~p} \textnormal{d}\Omega \ ... \nonumber \\ = \int_{\Omega}{~w^T\bar{~f}} \textnormal{d}\Omega + \oint_{\Gamma_N}{~w^T \bar{~h}} \textnormal{d}\Gamma \label{linwf1}\\
 \int_{\Omega}{~q^T\textbf{div}(~v)}\textnormal{d}\Omega  = 0 \\
 \nonumber
\end{align}

We can make the relevant substitutions for approximations as outlined above. This simplifies to a matrix system:

\begin{align}
\begin{array}{cc}
\sum{^{nel}_{e}}{\left(
\left[ \begin{array}{cc}
~c^T & ~q^T  \\
\end{array} \right]^e \left[ \begin{array}{cc}
\bar{~K}^{e} & ~G^{e}   \\
~G^{eT}  & ~0 \\
\end{array} \right] \left[ \begin{array}{c}
~\varphi  \\
~p  \\
\end{array} \right]^e  - \left[ \begin{array}{c}
~F^{e}  \\
~0  \\
\end{array} \right] \right)}
 &= 0\\
 \end{array}
\nonumber
\end{align} 

Where 
\begin{eqnarray}
\bar{~K}^{e} =& \int{~B^{e T} \bar{~D^{e}} ~B^{e}}\\
~G^{e}  =&\int{~L^{e T} \tilde{~N^{e}}}\\
~F^{e}  =& \int{~N^{e T}\bar{~f^{e}}}  + \oint{~N^{e T} \bar{~h^{e}}}
\end{eqnarray}

Which, through the arbitrary quality of the test functions and the nodal values independence on elemental summations, allows us to write a compact form:


\begin{equation}
\left[ \begin{array}{cc}
\bar{~K} & ~G  \\
~G^T & ~0 \\
\end{array} \right] \left[ \begin{array}{c}
~\varphi  \\
~p  \\
\end{array} \right] = \left[ \begin{array}{c}
~F \\
~0  \\
\end{array} \right]
\label{MM}
\end{equation}

\subsection{Solvability}
The system described above, with a null submatrix on the diagonal, governs velocity and pressure of our problem. We need to define under which conditions (\ref{MM}) can be safely solved. Provided that the velocity and pressure interpolations satisfy the $Babu\breve{s}ka-Brezzi$ compatibility condition\footnote{This condition ensures that the null space of $~G$ is zero \citep{stokeflow}}, they are uniquely defined and as a result the coefficient matrix (\ref{MM}) is non-singular and can be solved for. If one chooses an insufficient interpolation combination it can introduce spurious pressures or inadequately constrain the system.\\ 

\textbf{In the case of linear viscosity,} we would solve for $~v$ and $~p$ by partitioning the system and exploiting the positive-definite property of $\bar{K}$, allowing us to invert it:

\begin{eqnarray}
~\varphi  &=& \bar{~K}^{-1}(~F -~G~p ) \label{MMa} \\
~G^T ~\varphi &=& ~0 \label{MMb}\\
\nonumber
\end{eqnarray}

We substitute (\ref{MMa}) into (\ref{MMb}) which yields the solution for pressure:

\begin{equation}
(~G^T \bar{~K}^{-1}~G )~p  = ~G^T \bar{~K}^{-1} ~F \\
\label{MMc}
\end{equation}

Where the matrix $(~G^T \bar{~K}^{-1}~G )$ is symmetric and positive definite (granted that the null space of $~G$ is zero). Once the pressure values are solved, one can obtain the velocity result from (\ref{MMa}).

\textbf{In the non-linear case}, which this paper is focused on we would need to solve the system by the use of an iterative Newton method. To achieve this consider (\ref{linwf1}), rewritten as:

\begin{align}
\int_{\Omega}{ (\nabla_s ~w)^T \acute{~\sigma} (v)} \textnormal{d}\Omega + \int_{\Omega}{\textbf{div}(~w)^T ~p} \textnormal{d}\Omega \ ... \nonumber \\ = \int_{\Omega}{~w^T\bar{~f}} \textnormal{d}\Omega + \oint_{\Gamma_N}{~w^T \bar{~h}} \textnormal{d}\Gamma 
\label{eq:nl1}
\end{align}
 
Where $\acute{~\sigma}$ depends on $~v$, and therefore on the degrees of freedom $~\varphi$. Our system to be solved, (\ref{MM}) can be considered as:

\begin{equation}
\left[ \begin{array}{cc}
\bar{~K}(~\varphi) & ~G  \\
~G^T & ~0 \\
\end{array} \right] \left[ \begin{array}{c}
~\varphi  \\
~p  \\
\end{array} \right] = \left[ \begin{array}{c}
~F \\
~0  \\
\end{array} \right]
\label{MMnl}
\end{equation}

The tangential matrix, $K^{tan}$, used within in each iteration in the Newton scheme is then found from:
  
\begin{equation}
  K^{tan} = \frac{\partial}{\partial \varphi} [ \bar{~K}(~\varphi) ]\\
\label{eq:nl2}
\end{equation}

Which can be evaluated from the chain rule, and will be positive definite. For a specified load step, $\lambda_{k}$, the residual error of each iteration is found from:

\begin{equation}
R_{i} = \left[ \begin{array}{cc}
\bar{~K}(~\varphi_{i}) & ~G  \\
~G^T & ~0 \\
\end{array} \right] \left[ \begin{array}{c}
~\varphi_{i}  \\
~p_{i}  \\
\end{array} \right] - \lambda_{k} \left[ \begin{array}{c}
~F \\
~0  \\
\end{array} \right]
\label{MMnl2}
\end{equation}


Given initial starting vectors $~\varphi_{0}$ and $~p_{0}$, a loading step $\lambda_{k}$ we employ the Newton scheme to calculate in each iteration $K^{tan}_{i}$ and update the velocity and pressure vectors until the residual error is within a prescribed tolerance. The algorithm continues until a solution is reached for the system (\ref{MMnl}) with full loading.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Basal Terrain model}
To create a realistic bed topography we apply a numerical algorithm for creating a random terrain. As explained in Sec 3.1, we will only apply the bed topography to the centre of the bed in our model and preserve the smooth geometry leading to and emanating from our periodic boundary conditions (as depicted in fig \ref{fig:setup2}).\\

Given the domain of the terrain [0,L], the end points of the function $f(0)$ and $f(L)$, a perturbation maximum $h$ and a roughness parameter $r$ we calculate $f(mid)$ using a commonly used Random Midpoint Displacement (RMD) method, also know as Fractal Brownian Motion (FBM) Method \citep{Fournier}. Once $f(mid)$ is calculated we subdivide between known points and calculate new \textit{midpoint} values between points. \\

\begin{figure}[h]
	\centering
		\includegraphics[width=0.8\columnwidth]{pictures/3stepa.png}
	\caption{Basic results of RMD method. The number of known points after $n$ iterations is $2^n +1$.
}
	\label{fig:mid3}
\end{figure}

A key principle to the RMD is that with each further subdivision the random displacements should decrease by a certain factor which would depend on the prescribed regularization parameter $r$ which controls roughness. The most commonly used factor is given by $2^{-r}$, and $r$ is limited to the range $[0,1]$. As $r \rightarrow 1$, the factor $2^{-r} \rightarrow 1/2$ which corresponds to halving the potential displacement in each iteration, the largest possible ``smoothing'' that the RMD terrain generation allows.\\

The following result below in (fig.\ref{fig:Roughness}) tracks closely with our desired expectations.\\

\begin{figure}[h]
	\centering
		\includegraphics[width=1\columnwidth]{pictures/Roughness.png}
	\caption{Basic results of Midpoint displacement method for 8 subdivisions with varying roughness factor}
	\label{fig:Roughness}
\end{figure}

Once the satisfactory terrain was generated, it was imported as a set of displaced values into the model. I interpolated between them using a cubic spline function, to remove any kinks that may create discontinuities in the model. All finite element modelling is done in the \emph{COMSOL Multiphysics} environment, a commercial software package. The terrain generation modelling was created using \emph{MATLAB}.\\

 The following figure is the idealized reference geometry we will use to base our experiments on:

\begin{figure}[h]
	\centering
		\includegraphics[width=0.8\columnwidth]{pictures/mesh.png}
	\caption{Discretized 2D geometry with 756 nodes}
	\label{fig:2dmesh}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Basic Results}

This section will outline the solution to the Stokes flow problem, for a particular bed geometry created with our RMD code. If the solutions to the model seem feasible to the physical reality, we will use them in the following chapters as a controlled reference solution.\\

To induce velocities larger than $20myr^-1$ we choose the surface slope to be $\left| \alpha \right| = 0.2$ degrees.\\

To start, we see that the basic magnitude for velocity seems feasible. Firstly the range of our velocity magnitude falls in line with generally accepted levels outlined in `\textit{Principles of Glacier mechanics}' \citep{pgm}. We are particularly interested in surface velocity, which we will use in the forthcoming sections as a means of comparing optimization successes.\\
  
 \begin{figure}[H]
  \centering
  	\subfloat[Magnitude of velocity]{\includegraphics[width=0.45\columnwidth]{pictures/velmag.png}}                
  	   \subfloat[Surface velocity]{{ }\includegraphics[width=0.45\columnwidth]{pictures/velmagS.png}}\\
  \caption{Reference configuration Velocity Magnitude}
  \label{fig:referenceV}
\end{figure}


We introduce a scaling measure on our basal traction, $\beta^{2}_{opt}$, which will be the variable we will optimize in the later experiments. This changes our boundary condition specifying traction (\ref{eq:bedslide}) to: 

\begin{equation}
~\sigma \cdot ~n = \beta^{2}_{opt} \beta^2 ~v 
\label{eq:bedslide2}
\end{equation}

Where $\beta^2$ has been set to $1000 Pa.yr.m^{-1}$ and $\beta^{2}_{opt}$ is trivially equal to 1 over the entire bed for the reference configuration. In the following chapter of optimization, we make $\beta^{2}_{opt}$ the focus of our minimization routine and in this way directly scale the traction of our boundary condition (\ref{eq:bedslide2}).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Optimization}

At this stage we return to our posed problem with the lack of accurate bed data. We take our \emph{Reference} model that was described in the preceding chapter and present it as a theoretical perfect solution. This model has been solved for with our RMD generated terrain using 33 points ($2^5$ divisions). We store the surface velocity as a means of comparison. We deliberately introduce secondary models, with incorrect bed geometries. We wish to adapt variables within the secondary models so that they match the behaviour of our correct, reference model.\\

We take the squared difference between the surface velocities of our reference and secondary models as an error in the secondary models that we wish to remove. To achieve this we employ a numerical optimization algorithm, varying the scaling value for the traction field, $\beta^{2}_{opt}$, at the ice-rock interface to minimize the error.\\

Due to the manner in which we created our bed geometry, with smooth linear bed geometries preceding and leading to the vertical boundaries. We expect our defined error to be concentrated over the centre of our domain as illustrated in the figure below:\\

\begin{figure}[H]
	\centering
		\includegraphics[width=.6\columnwidth]{pictures/objective.png}
	\caption{The figure above illustrates the error in surface velocities between a reference geometry solution and a perturbed bed geometry solution}
	\label{fig:objective}
\end{figure}

Mathematically we can express the optimization problem in terms of:\\
\begin{eqnarray}
\nonumber	\text{minimize  } & \ & F(~x) \\
\text{subject to  } & \ & c_{i}(~x) = 0 \text{ \ where $i = 1,2,..m$}
							\label{min}
\end{eqnarray}

The objective function to be minimized, $F$, is defined by all input parameters, $~x$, from the model. The constraints on the parameters, $c_{i}(~x)$, for $m$ constraints determine variables (and the extent to which) the optimization algorithm can manipulate. For our specific case we set the $\beta^{2}_{opt}$ field as our sole control variable for the optimization algorithm, strictly constraining all other input parameters that determine the objective function, the error between the surface velocities.\\

\subsection{Tikhonov regularization}
It was found that initial solutions to the posed optimization problem favoured high frequency solutions, demonstrated in (fig: \ref{fig:sampleBopt}) , which in turn model rapidly changing ice flow behavior on the bed. Supporting the hypothesis that errors exist in our basal topogrophy (introduced be an error in approximating geometry or traction) and that these errors in parameters affect ice flow, it follows that smooth solutions would more accurately represent the true behavior of the ice flow.\\

 \begin{figure}[H]
  \centering
 \subfloat[$\beta^{2}_{opt}$]{\includegraphics[width=0.45\columnwidth]{pictures/BetaOpt.png}}                
 \subfloat[Velocity profile of ice sheet]{{ }\includegraphics[width=0.45\columnwidth]{pictures/Velmagaft.png}}\\
  \caption{Scaling factor $\beta^{2}_{opt}$, solved for in the optimization routine and influence on ice flow}
  \label{fig:sampleBopt}
\end{figure}

To facilitate realistic, smooth solutions we introduce a Tikhonov regularization parameter \citep{Tink}, such that (\ref{min}) becomes:\\

\begin{eqnarray}
\nonumber	\text{minimize  } & \ & \lbrace F(~x) +\lambda \Vert ~L x \Vert ^2 \rbrace \\
\text{subject to  } & \ & c_{i}(~x) = 0 \text{ \ where $i = 1,2,..m$}
							\label{min2}
\end{eqnarray}

Where $\lambda $ is a scaling parameter and $~L$ is the Tikhonov regularization matrix. For our specific problem $~L$ is set to be a first order horizontal gradient operator. The addition of the term $\lambda \Vert ~L x \Vert ^2$ to the minimization routine penalizes solutions with high frequency and favours those that are smoother.\\

\emph{COMSOL} uses the SNOPT optimization algorithm; a sequential quadratic programming (SQP) method developed by \textit{P.E Gill, W Murray and M.A Saunders}. The method solves nonlinear problems by iterating through sequences of quadratic subproblems (QP). In our case, we set SNOPT to use the \textit{Quasi Newton} subroutine to solve our minimization problem and set the algorithm to run until the error is reduced to be smaller than $1 \times 10^{-4}$. A detailed explanation of the SNOPT routine can be found in the paper \textit{SNOPT: An SQP algorithm for large-Scale Constrained Optimization}\citep{snopt}.\\

A thorough guide to optimization can be found in the books \textit{Numerical optimization} \citep{nocedal} and \textit{Practical Optimization} \citep{opt1}.\\
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Further Results \& Analysis}

In this chapter we introduce perturbed bed geometries to be tested by the optimization routine. In each case, the basal friction will be manipulated, via the optimized $\beta^{2}_{opt}$ scalar field, to force the model's surface velocity to conform with our reference model.\\

We will focus the experiments to the perturbation of a single point along the bed geometry of increasing size. Each experiment will increase the error in the bed geometry and record the resulting \textbf{Flux divergence/particle tracing} which will be compared to the original ``true'' reference case. \\


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\conclusions[Conclusions and Future directions]
TBC
%Taking our top three performing cases and plotting them against each other confirms the findings in the previous sections. Interestingly, we can see that the error is reduced uniformly in all cases at the centre of the surface, corresponding to the lowest values of surface velocity. Additionally, there are distinct points where all three solutions have the same error value. One must assume that these are points isolated and weighted by the optimization algorithm. 
%
%\begin{figure} [H]
%	\centering
%		\includegraphics[width=.8\columnwidth]{pictures/thebest/Error.png}
%	\caption{Error in the top performing cases, aligned to the Reference Velocity}
%	\label{fig:Error}
%\end{figure}
%
%The accuracy of these solutions, measured point wise, has an accuracy greater than $99.5\%$. While the flat bed scenrio may not be as effective as the above results, once could allow for it in analysis which did not require a high evel of precision.\\
%
%The largest implication is in the use for models that have data for bed elevation but one suspects inaccuracies within the data... Perhaps the data resolution has allowed for a large dip or rise to have been smoothed over too dramatically, or perhaps the data itself is mildly inaccurate do to the nature of the radar reading equipment. In these cases the optimization routine employed here would benefit the models and allow for more accurate results, with the manipulation of the bed data itself.\\
%
%\subsection{Future direction}
%The next step in this research is to investigate how significantly the ice thickness evolution is affected by the optimization. Ice thickness, $H$, is described by the flux divergence equation:
%
%\begin{equation}
%\frac{\partial H}{\partial t} = ~M - \textbf{div}\left( \bar{U} H \right)
%\label{eq:height}
%\end{equation}
%
%Where $H$ is the ice thickness, $~M$ is the net accumulation rate and $\bar{~U}$ is the depth averaged velocity. An essential question of glacier study is whether or not a glacier range is growing or shrinking. An accurate solution to the ice thickness evolution is an important measure. Whether or not our solution to the optimization problem could be fully accepted would be dependant on how radically the ice thickness evolution changes.\\



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliographystyle{cell}
\bibliography{mybibliography}{}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\end{document}

